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The process of verification, or determining the order of accuracy of computational codes, 
can be problematic when working with large, legacy computational methods that have been 
used extensively in industry or government. Verification does not ensure that the computer 
program is producing a physically correct solution, it ensures merely that the observed 
order of accuracy of solutions are the same as the theoretical order of accuracy. The 
Method of Manufactured Solutions (MMS) is one of several ways for determining the order 
of accuracy. MMS is used to verify a series of computer codes progressing in sophistication 
from “textbook” to “real life” applications. The degree of numerical precision in the 
computations considerably influenced the range of mesh density to achieve the theoretical 
order of accuracy even for 1-D problems. The choice of manufactured solutions and mesh 
form shifted the observed order in specific areas but not in general. Solution residual 
(iterative) convergence was not always achieved for 2-D Euler manufactured solutions. 
L 2 , norm convergence differed variable to variable therefore an observed order of accuracy 
could not be determined conclusively in all cases, the cause of which is currently under 
investigation. 


N omenclature 


S Source term 

C Coefficient in manufactured solution 

E Energy 

L,L x ,L y Reference length 

N Number of mesh nodes between computational boundaries 

p Pressure 

Re 1/v 

t Time 

u, v Cartesian velocity components 

x, y Cartesian coordinates in physical space 

L 2 ,norm vector norm 

T.E.i (Ax+ — Ax-) 4>xxx ) First term of truncation error, Eq. 12 

1 f ~\~ ^^X^ \ 

T.E.a j 2 ( Ax++Ax~ ) ( t ) xxxx i Second term of truncation error, Eq. 12 
Subscripts 

0 Initial mesh spacing at x=0 

MS Manufactured solution 

1 Node index 

Symbols 

7 gas constant, 1.4 

Ax Inter-node spacing 
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Ax+ (x i+ 1 - Xi) 

Ax _ (xi - i) 

Ax n Equal interval spacing, 1 /N + 1 

v Kinematic viscosity 

f 1 Computational domain 

dtt Boundary of computational domain 

0 Arbitrary function in physical space 

p Density 


I. Introduction 


A broad definition of the numerical verification of a program written to solve partial differential equations 
(PDE’s) is succinctly stated by Knupp & Salari 1 

“The process by which one demonstrates that a PDE code correctly solves its governing equa- 
tions.” 


Consider the Laplace equation, V 2 0 = 0, within some boundary Q with a set of conditions on the 
boundary, 0 = g. For this discussion we consider solutions using finite difference techniques rather than 

any of the many analytic techniques that exist for this equation. Discretized equations are approximations of 
derivatives arrived at typically through Taylor series expansions of, in this example, the dependent variable 
0 at some distance Ax from x. One discrete expression for the first derivative of 0, with order property Ax , 
can be written as, 


d 0 ^ (j){x + Ax) — (j){x ) 
dx Ax 


O ( Ax ) 


For Laplace’s equation, an expression for the second derivative could be written as 


(1) 


d 2 0 0(x + Ax) — 20(x) + 4>(x — Ax) 

dx 2 Ax 2 

That is, as Ax is decreased, the discrete expression (0(x + Ax) — 20(x) + 0(x — Ax)/ Ax 2 should approach the 
exact expression d 2 0/dx 2 at the rate of Ax 2 . It is important therefore to verify that a program, written for 
example using Eq. 2 in the discretization of the second derivative, exhibit the correct second order property 
behavior. 

The method of manufactured solutions is one of several ways to verify the order properties of numerical 
methods. Discussions of verification methods can be found in various books and reports, for example, by 
Roache 2 , Salari & Knupp 3 , Oberkampf & Trucano 4 and Roy et al. 5 
Briefly, consider the general form, 


+ 0(Ax 2 ) 


( 2 ) 


V<j)-S = 0 


( 3 ) 


where V and S are, respectively, a differential operator on 0 and some source term. When using the method 
of manufactured solutions for verification, another source term, iS MS , is added to Eq. 3 as a forcing function 


V(f> - S = 5 ms , where 5 MS 


V<j)-S 


4 * — 0MS 


( 4 ) 


The function 5 MS is a function consisting of the same differential operator as the modeled equation applied to 
the predetermined solution function, 0 ms- The difference between the iteratively converged discrete solution 
of Eq. 4 and the predetermined solution function is the discretization error and will be a function of the 
mesh spacing and the truncation error of the discretization scheme. The sum of the differences between the 
discrete solution and the predetermined solution over the solution space results in an integrated value for 
the discretization error L 2 , n orm a 


L2,] 


(0 solution 0MS ) 

N 


( 5 ) 


This process is repeated for 3 or more mesh densities. The observed order property of the solution method 
is obtained from the change of the discretization error with mesh size. 


a http: //mathworld. wolfram.com/L2-Norm.html 
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II. 1-D Poisson Equation 


The ubiquitous Poisson equation subject to Dirichlet boundary conditions, Eq. 6, serves as the initial 
platform for demonstrating text book verification, both successes and pitfalls, using manufactured solutions. 
The function f satisfies Eq. 6 throughout the interior of fl subject to the physical source term v. The 
Dirichlet condition f = g is set along the boundary dfl. 

V 2 f — v = 5ms, in D, f = g, on dfl. (6) 

The manufactured solution source term is defined in Eq. 7, where f M s and u MS are specified functions. 


5ms — V Ims ^ms 


( 7 ) 


Now consider the domain fi to be a single dimension partitioned in to N segments with the boundaries at 
x 0 and x N+1 such that x 0 < Xi < x N+1 , 0 < i < N + 1. 

II. A. Three point discretization in physical space 

If / is a function of in physical space x, the difference of slopes at the half nodes z ± 1/2 can be used to 
calculate the second derivative of /. 


(8) 


4*x | 

i+l/2 

2 

( 02+1 02 02 02—1 \ 

^i+l/2 %i- 

-i /2 Ax+ + Ax- 

iy Ax. Ax- ) 


using the averaging expressions 


U+l/2 


4>i + 1 

Ax . |_ 


* 1 / 1 Ax- 


(9) 


where (f>i = </>+), Ax+ = Xi+i — Xi , and Ax- = Xi — i . The Taylor series expansions of <f) (Eq. 10) is 

substituted in to Eq. 8 


^z±l — ((ii | T l ^(j) X x\^Ax‘^. i ^fpxxxl^Ax^. “b ^^4>xxxx\ ^Ax± 

190 ^ xxxxx T ~j2Q ( t >x xxxxx\i Ax± i 0 (Ax±) 

resulting in the expression for f xx (Eq. 11). 

f xx |j = 4*xx T 


( 10 ) 


( 11 ) 


where truncation error terms of f x 


IT.E. 


are 


, I 1 . . . . | If Ax 3 + Ax 3 * _ 

fxx |t.E. - 3^*+ Ax -)<t>**x\i+ 12 [ Ax+ + Ax- 


1 / Ax 5 + + Ax 5 _ 


360 \ Ax+ + Ax- 


1 / Ax\ — Axf 
* + 60 \ Ax + + Ax- _ 

LXXXXX l + 0(Ax 5 ± ) 


( 12 ) 


If the mesh is uniformly distributed so that Ax = Ax + = Ax_, then Eqs. 11 and 12 reduce to the well 
known forms 


fxxlj — 


4>i + 1 — 2(/>j + <fii- 1 
Ax 1 


and 


fx 


IT.E. 


12 


^ Ax ^qq^xxxxxx\ ^A x + O (Ax_$f) 


(13) 


(14) 


Equation 11 is formally only first order accurate in a non- uniformly spaced mesh since the leading term 
of Eq. 12 goes as Ax. In a uniformly spaced mesh, the leading term of Eq. 14 varies as Zia; 2 so that the 
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difference equation Eq. 13 is formally second order accurate. Though, as noted in Hoffman 6 , if the quantity 
(Ax+ — Ax-) has the characteristic of decreasing by a factor of 4 for a halving of the mesh spacing, then 
Eq. 11 is second order accurate as well since the leading term in Eq. 12 would then decrease as 0 (Z\a; 2 ). 

Also discussed in Hoffman, 6 but without detailed derivations, were expressions for the central difference 
expression for the second derivative in transformed space. The derivation is similar to the one in physical 
space just discussed and is detailed in Appendix A, 

II. B. Forms of Manufactured Solutions 

Several guidelines have been published concerning the form of the manufactured solutions to be used in 
verification of computational methods, see Roache 2 or Knupp & Salari. 1 As far as the assessment of a 
methods order of accuracy is concerned, the most important feature for any manufactured solution is the 
truncation error should be non-zero. 


II.B.l. Exponential and Trigonometric Functions 

Exponential and trigonometric functions are often used in manufactured solutions. These functions are 
smooth, continuous with an infinite number of smooth and continuous higher order derivatives. 

A simple example would be single exponential function 

jm f 

/ms = e-* - = (-ire-* (15) 

as do trigonometric functions such as sin or cos, and hyperbolic functions such as sink and tank. For 
example , 


n • /2ttx 
= ( sin 


/ms = c sin 


V L 


d™f MS = \C(-l)^m m cosm : modd 

dx m | C(-ir (x) m sin(^f) : m even 


(16) 


Using a sin function as the manufactured solution, Eq. 16, the MS source term is 


5ms — — C 


2-7T 


. ( 2nx \ 
in — 

V L ) 


sin 


%s 


(17) 


II. B. 2. Polynomial Series 


It is also possible to use a polynomial in x to some power n as a manufactured solution, for example Eq. 19. 


/ms = x n 


[/msU = n(n - 1) ; 


v MS = v(x) 

The source term for the manufactured solution is then 


<Sv 


= n(n - 1) x” - v MS 


(18) 

(19) 

(20) 


As mentioned previously, the manufactured solutions should be chosen such that the truncation error is 
non-zero for a finite Ax in the modeled equations. A polynomial manufactured solution has to be of power 
n = 4 or higher for a uniformly spaced mesh since the leading term of the truncation error is 1/12 Ax 2 f X xxx- 
If n < 4 the manufactured solution will be an exact solution of the discretized equations and not display 
any order property with mesh refinement. There is a potential benefit from this result and is discussed in a 
subsequent section. 


II. C. 1-D Verification Case Studies 

Changes ins L 2 jn0 rm of the observed order and truncation error with the form of manufactured solution, 
mesh distribution, stencil size and numerical precision using the 1-D Poisson equation is discussed in the 
following sections. The parameter variations are listed in Table 1; though not all the cases will be specifically 
discussed. 
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Table 1. List of Case Studies 


Case 

/ms 

Stencil 

Mesh Form 

1 

sin(x ) 

3 pt. 

Uniform 

2 

sin{x ) 

3 pt. 

Geometric series 

3 

x 5 6 

3 pt. 

Uniform 

4 

x 3 

3 pt. 

Uniform 

5 

X 3 

3 pt. 

Geometric series 

6 

iV 7^ V MS 

3 pt. 

Uniform 

7 

X 3 ,V ± V MS 

3 pt. 

Uniform 

8 

X 3 

3 pt. 

Elliptic DE * 

9 

X 6 

3 pt. 

Elliptic DE 

10 

sin(x) 

5 pt. 

Uniform 


Differential Equation 


II.C.l. Study 1 - A Textbook Solution 


Table 2. Study 1 Parameters 


Case 

Symbol 

/ms 

%S 

Mesh 

real*4 

O 

sin( 2ttx) 

V 

Uniform 

real*8 

□ 

sin{ 2irx) 

V 

Uniform 


Consider a uniformly space computational 
mesh, Ax with Dirichlet conditions applied 
at the boundary points, <p\ x=0 = 4>\ X=L = 
0, and a manufactured solution of the form 
sin (2-kx/L). Single (32-bit) and double (64- 
bit) precision calculations were performed us- 
ing a direct solver on the second order accu- 
rate 3 point stencil equations in physical space 
shown in Eq. 21. The precision of the pro- 
gram variables were explicitly defined using the 
real*4 or real*8 FORTRAN statement calls 
for the single and double precision calculations 
respectively. The source code was compiled 
with Intel F90 FORTRAN and executed on the 
1.4 GHz Opteron CPU. 




(a) I j 2,norm vs. 1/N (b) Observed order vs. 1/N 

Figure 1. 1-D Poisson equation, 3-point stencil, 2 nd order scheme, 
O? 32-bit precision ; □, 64-bit precision 


4>i + 1 — 2 (f)i + 4>i-l — 

— 4-7T 2 sin (2-KXi) Ax 2 , 't , Ms| i = 'U 


( 21 ) 


The variation of the discretization error, I^norm, with mesh spacing and the observed order of accuracy 
is shown in figure 1. The number of computational nodes ranged from 5 to 540,217 corresponding to a 
node-to-node spacing from approximately 1.7 xlO -1 to 1.9xl0 -6 based on a unit interval. The range also 
corresponds to a mesh ratio of r=1.5 for 30 mesh levels starting at N=5. 

The design numerical order using Single precision arithmetic was limited in the range of mesh spacing 
due to accumulative round-off error. The use of double-precision arithmetic allowed for another 2.5 orders 
of magnitude reduction in mesh spacing and still retain textbook order of accuracy, figure 1(b). 


5 of 18 


American Institute of Aeronautics and Astronautics Paper 2006-7128 


II. C. 2. Study 2 - Difference with geometric series mesh expansion , Case 1 vs. Case 2 


Table 3. Study 2 - Case 1 vs. Case 2 


Case 

Symbol 

/ms 

%S 

Mesh 

T.E.i 

T.E .2 

1 

O 

sin( 2 ttx) 

V 

Uniform 

0 

-b (27r) 4 Z\a; 2 S'm(27ra;) 

2 

□,o 

sin( 2irx) 

V 

Geometric 

series 

— |(27t) 3 ( Ax+ — Ax-) cos{2t:x) 

12 W ( Ax + + +Ax-) sin{2nx) 


The mesh spacing for case 1 was 
Ax = 1 /(TV + 1) between nodes. 
For case 2 the spacing interval in- 
creased away from x = 0 with the 
first interval Ax o = 0.0001/(iV+ 1) 
for each of the meshes. A Newton 
iteration method was used to de- 
termine the correct expansion rate, 
r, to fit the N nodes between the 
boundaries 0 < x < 1 as defined in 
Eq. 22. For the cases studies here 
L = 1. 

/ N \ 




1/N 


L = Ax i 



( 22 ) 


(a) t j 2 ,norm vs. 1/N , O, Uniform mesh ; (b) Observed Order vs. 1/N , O, Uni- 

□ , Geometric series ; , 0(Ax 2 ). form mesh ; □, Geometric series ; , 

O {Ax 2 ) . 


The individual node positions were 
calculated using Eq. 23 

Xi+i = Xi + Axor 1 (23) 

Table 10 in Appendix B details the 
total number of nodes, initial spac- 
ing and geometric expansion ratio 
for each mesh run. A direct solver 
was used to calculate solutions with 
Dirichlet boundary conditions for 
all 30 mesh density levels. 

The variation of the discretiza- 
tion error, I^normi with mesh spac- 
ing and the observed order of ac- 
curacy for the two mesh distribu- 
tion are compared in figures 2(a) 
and 2(b). Figure 2(c) is a plot of 
the Li 2 ,t.e. of the non-zero trunca- 
tion error terms for the 2 cases and 
for the N = 163 geometric series, the first two terms of truncation error mesh are shown in figure 2(d). 

The uniform and geometric series mesh for Kt 100 < N <« 10, 000, both displayed with second order 
accurate behavior using a sinusoidal manufactured solution, comparing the O and □ symbols in figures 2(a) 
and 2(b). The asymptotic range was roughly the same between the two mesh types-the uniformly spaced 
mesh approached order 2 in a more continuous fashion than the non-uniform meshes and departed from the 
design order due to numerical round off errors at a lower value of N. The different trends of the L 2 with 
increasing N is due to significantly higher truncation error of the non-uniform grid for meshes with less than 
70 nodes. For N=73 the geometric expansion factor was approximately 1.2, see Table 10. 

The L 2 , norm of T.E .2 were almost identical between the two cases, compare O and O symbols in figure 2(c). 
The first and second truncation error terms for case 2 are similar in magnitude, but 7 r radians out of phase, 



(c) L 2 y T.E vs. 1/N, O, Case 1 - T.E .2 ; 
□ , Case 2 - T.E.i ; O, Case 2 - T.E .2 ; 
, O {Ax 2 ) 

Figure 2. 1-D Poisson equation, 

non-uniform mesh . 



(d) fxx, T.E.i vs. x, N = 163, ; -, 

i X x discrete solution ; , T.E.i ; , 

T.E.2 

sinusoidal MS, /ms — sin(x), uniform vs. 
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compare with — • — ■ in figure 2(d). The two truncation terms predominate in different regions of the 

solution space but overall the observed order is unaffected because both terms behave as Ax 2 in this case. , 
compare □ and O symbols in figure 2(c). 

II. C. 3. Study 3 - Polynomial manufactured solutions on a uniform mesh, Case 3 vs. Case 4 


Table 4. Study 3 - Case 3 vs. Case 4 


Case 

Symbol 

/ms 

%s 

Mesh 

T.E.i 

T.E.2 

3 

O 

X 6 

V 

Uniform 

0 

30Z\x 2 x 2 

4 

□ 

X 3 

V 

Uniform 

0 

0 


On the uniformly spaced mesh, 
the difference in the L 2 )nor m and 
the observed order between an or- 
der 6 polynomial and an order 3 
polynomial manufactured solution 
is shown in figure 3. In the former 
case the MS source term will be an 
order 4 polynomial, while the latter 
choice the MS a linear function of 
x. 

The first term of the truncation 
error, T.E.i, is zero for both cases 
due to the uniformly spaced mesh 
and only case 3 has a non-zero sec- 
ond truncation term, T.E. 2 , since 
the fourth derivative of x 3 is zero. 

Therefore the correct order proper- 
ties are seen only for the /ms = x 6 

manufactured solution, see the O symbol in figures 3(a) and 3(b). 

The function /ms = x 3 as a manufactured solution is an exact solution of the discrete equations, hence 
the □ symbol in figure 3(a) is representative of only accumulated round off error for L 2 jn0 rm and as a result 
in figure 3(b), does not display an appropriate observed order. 




(a) L 2 ,n 


l/N 


Figure 3. 1-D Poisson equation. 

, 0(Ax 2 ) . 


(b) Observed Order vs. l/N . 
uniform mesh, O, /ms = x 6 ; □, /ms = x 3 


II. C. 4- Study 4 ~ Polynomial manufactured solutions on a geometric series mesh, Case 4 vs. Case 5 


Table 5. Study 4 - Case 4 vs. Case 5 


Case 

Symbol 

/ms 

V MS 

Mesh 

T.E.i 

T.E.2 

4 

O 

X 3 

V 

Uniform 

0 

0 

5 

□ 

X 3 

V 

Geometric series 

2 (Ax+ — Ax-) 

0 


Case 5 used the same manufactured solution as did Case 4 but was solved on a non-uniformly spaced 
(geometric series ) mesh. An interesting result of this is that the first truncation error term, T.E.i, is non- 
zero and so an observed order should now be seen using a polynomial manufactured solution as low as order 
3. 

The L 2j norm and observed order for cases 4 and 5 are compared in figures 4(a) and 4(b). The L 2 of the 
first truncation error term, T.E.i, for case 5 is plotted in figure 4(c). 

Similar to what was discussed previously in section II. C. 2, the quantity ( Ax+ — Ax _) is non-zero and has 
the appropriate behavior to give the term T.E.i a second order convergence property with mesh refinement. 
Similar behavior seen in □ symbols of figure 4. 
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II. C. 5. Study 5 - Effect of physical source term error, Case 3 vs. Case 6 



(a) I j 2 1 norm vs. 1/N . (b) Observed Order vs. 1/N . (c) L 2 ,t.e.i vs. 1/N . 

Figure 4. 1-D Poisson equation, polynomial MS, /ms = x 3 , O, Uniform mesh ; □, Geometric series mesh ; 

, O (Ax 2 ) . 


Table 6. Study 5 - Case 3 vs. Case 6 


Case 

Symbol 

/ms 

%s 

Mesh 

T.E.i 

T.E.2 

3 

O 

X 6 

e~ x 

Uniform 

0 

30Z\ie 2 :e 2 

6 

□ 

X 6 

e -1.00001x 

Uniform 

0 

30Z\a; 2 ir 2 

7 

o 

X 3 

e -1.00001a; 

Uniform 

0 

0 


This section dis- 
cusses introducing 
an error in the 
physical source term 
function v of Eq. 6. 

Previously, the man- 
ufactured solution 
i>ms was identically 
v, so the source 
term was effectively 
“zeroed out” of the 

discretized equations, (a) L 2 , n orm vs. 1/N, / M s = 
The case 6 and case °> *ms = v , n, »ms ¥= v ; 




(b) Observed Order vs. 1/N, 
/ms = x 6 , o, U M S = V , D , 
»ms ^ V ; > o (Ax 2 ). 


(c) L 2 ,norm VS. 1/N, f M S = X 3 
O, »MS V ; , 0{Ax 2 ). 


O (Ax 2 ). 

Figure 5. 1-D Poisson equation, uniform mesh, polynomial MS, vms = v vs. i»ms 7^ v 


3 manufactured so- 
lution, /ms, and 
meshes in this sec- 
tion were the same. The manufactured physical source term -ums was modified to mimic an error in v, shown 
in Table 6. The L 2 , n orm and observed order for cases 3 and 6 are compared in figures 5(a) and 5(b). 

It is interesting to note that if the number of nodes in the order of accuracy study had been not exceeded 
1000, b the departure from second order behavior for case 6 would not have been observed. For comparison, 
case 6 was repeated with an order 3 polynomial for the manufactured solution instead of the order 6 poly- 
nomial. It is known from the discussion in section II. C. 3 that this combination of MS and mesh topology 
should produce I^norm values close to floating point zero plus accumulated round off error. The significantly 
higher and constant L^norm shown in figure 5(c) is an indicator of the physical source term error. 

Therefore the discontinuity in the I^norm distribution around 1/N ~ 10~ 3 for case 6, □ symbols in 
figure 5(a), is due to the physical model error. The same effect is observed in figure 5(b) with the shift from 
the correct observed order 2 down to 0 for the □. 


*’Tlie threshold of 1000 nodes is somewhat arbitrary in this case. The difference in the level of the L 2 , n 
equation compared with the I^norm induced by the physical source term error would determine this. 


of the discretized 
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II. C. 6. Study 6 - Polynomial manufactured solutions on a elliptic DE generated mesh , Case 3 vs. Case 9 


Table 7. Study 6 - Case 3 vs. Case 9 


Case 

Symbol 

/ms 

V MS 

Mesh 

T.E.i 

T.E.2 

3 

O 

X 6 

V 

Uniform 

0 

30Z\x 2 x 2 

9 

□ 

X 6 

V 

Elliptic DE 

40 ( Ax + — Ax-)x 3 




Iff 6 Iff 5 Iff 4 Iff 3 Iff 2 Iff 1 10° 
1/N 



The non-uniform 
mesh distribution 
results shown in 
sections II. C. 2 and 
II.C.4 had the cor- 
rect observed order 
due to the second 
order behavior of 
{Ax+ — Ax-) with 
mesh refinement. The 
non-uniform mesh 
for Case 9 was 
generated using an 
elliptic differential 
equation, Eq. 24, 
instead of the ge- 
ometric relation used 

in the previous sections. The coefficients of the function P were chosen so as to cluster the nodes closer to 
x = 0. 


l/N 



(a) L 2 ,norm vs. 1/N , O, Uniform (b) Observed Order vs. 1/N , O, (c) \j 2 ,t.e vs. 1/N, O, Case 3 - 


mesh , □, Elliptic DE mesh ; 
Q(Ax 2 ) . 


Uniform mesh , □, Elliptic DE T.E .2 ; n, Case 9 - T.E.i ; O, 


mesh ; , O (Ax 2 ) . 


Case 9 - T.E .2 


O {Ax 2 ) 


Figure 6. 1-D Poisson equation, polynomial MS, /ms — x& ? uniform vs. elliptic DE generated 

mesh . 


d 2 x 
d £ ^ 


„ dx 

P-77 = 0 

df 


The finite difference form used to generate the mesh and the clustering function was 


R 


R 


1 - — Xi-i - 2xi + I 1 + — x i+ 1 =0, Pi = - 


10 


o -i/(10(JV+l)) 


N + l 


(24) 


(25) 


Table 7 summarizes the parameters of this comparison study and Table 11 in Appendix B lists the initial 
mesh spacing with total node count. The L 2 . n orm and observed order for cases 3 and 9 are compared in 
figures 6(a) and 6(b). The L 2 of the truncation error terms are shown in figure 6(c). 

The coefficients used in the clustering function P were not optimized for any particular mesh. As can be 
seen in figures 6(a) and 6(b), □ symbols, meshes with less than 50 nodes were not in the asymptotic range 
as the L 2 ,norm distribution was nearly flat and the observed order was 1.9 or less. The clustering function 
was excessive for the more sparse meshes and could have been adjusted, though as far as this case study was 
concerned, wasn’t really required as the design order of accuracy was still observed, just within a smaller 
range of mesh densities. 

The relevant aspect of the truncation error terms for case 9, plotted in figure 6(c), is that both T.E.i and 
T.E .2 have second order characteristics with increasing mesh density for N/50. The elliptically generated 
mesh extended the range of asymptotic convergence to higher values of N than with the uniformly spaced 
mesh. This result is likely academic as the particular node count and spacing would be unrealistic for most 
3-D meshes. 
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II. C. 7. Study 7 - A five-point stencil for the second derivative - Case 10 


Table 8. Study 7 Parameters 


Case 

Symbol 

/ms 

%S 

Mesh 

real*4 

O 

sin(2irx) 

V 

Uniform 

real*8 

□ 

sin(2irx) 

V 

Uniform 


The final case study for the 1-D Poisson problem investigated the order 
characteristics of a higher order accurate stencil. Similar to section II.C.l, 
single and double precision arithmetic computations are compared over 
a range of mesh densities. The L 2 . n orm and observed order are shown in 
figures 7(a) and 7(b). 

A penta-diagonal solver was used to resolve the discrete set of equa- 
tions from applying a central 5 point finite difference stencil , Eq. 26, to 
a uniformly spaced physical grid. 

- (fi -2 + 16^_i - 30(fi + lG4>i + i - (fi+2 = RHS| . (26) 

Explicit boundary conditions were imposed at the end points modifying 
the first two and last two discretized equations were modified as shown in 
Eq. 27. 


RHS\ 1 


125 | j Ax N + <^ms | _i 16</> ms | Q 

rhs\ 2 


12 S\ 2 Ax 2 N + ^ MS 0 

RHS\ 

\ l 


12 S\.Ax% 

RHS \n-i 


125 \ N -! Ax2 n + />ms |jv + i 

rhs\ n , 

= 

• \n^ x N - 160ms 1^1 + </ms 1 jv_|_2 ' 


It was surprising to observe the rather rapid departure from the de- 
sign order of 4 with the use of single precision arithmetic, O symbol in 
figures 7(a) and 7(b). 

Double-precision arithmetic allowed for almost a 4 order reduction in 
mesh spacing and still retain the textbook order of accuracy, □ symbol 
in figure 7(b). So by comparison, the higher order stencil had a smaller 
range of applicability while retaining the design order of accuracy of the 
discretization scheme. 



(a) L 2 ,norm vs. 1 /N . 4 th order scheme. 



1/N 


(b) Observed order vs. 1/N , 4 th order 
scheme. 

Figure 7. 1-D Poisson equation, 5 

point stencil, effect of machine pre- 
cision, Oj single precision ; □, dou- 
ble precision . 


III. 2-D Burger’s Equation 


Advancing in complexity, the second “textbook” problem was the 
advective-diffusive Burger’s equation, where u and v are velocity com- 
ponents, and Re is the inverse of the kinematic viscosity. 


du du If d 2 u d 2 u\ 

U dx V dy Re \ dx 2 dy 2 ) 

dv dv if d 2 v d 2 v\ 

U dx V dy Re \ dx 2 dy 2 ) 


(28) 

(29) 


The manufactured solutions and source terms had the form: 


u MS = u 0 + Sill 


. ( 2ttx\ ( 27 xy 

m I I -i- cos I 


V L x J 


\ Ly 


, V MS = v 0 + cos 


f 27Tx\ 

lM 



(30) 
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, MS 


(31) 


5«,MS — 
S V .MS 


du 

du 

71 _ 

_2_I 

( d 2 u 

d 2 u 

dx 

dy 

Re 1 

\Jh? 

+ dy 2 

dv 

dv 

A m 

1 , 

f d 2 v 

d 2 v' 

dx 

+ — 
dy 

Re 1 

\dx 2 

+ dy 2 . 


' 


»=«MS,»=»MS 


(32) 


where uq = vq = 0.5, L x = L y = 7t/( 6A), A = 25, and Re = 50. 

An 2-factor approximately factored scheme with 2 nd order, central- 
differencing was used to discretize the equations with pseudo-time stepping to 
obtain iterative convergence. The system of equations were solved using the 
cyclic Thomas algorithm along each direction. The computational domain 
was 2-dimensional with an even spaced mesh in the vertical and horizontal 
directions. Dirichlet boundary conditions were imposed on the four edges of 
the domain, u = v = 0. This computer code approached the appropriate 
order property as the inverse mesh density, 1/N = 1 /{N x xJV 9 ), approached 
0.0001 denoted as the “asymptotic range” in figure 8. The data below that 
of 2 is an indication of too coarse a mesh for the numerical scheme to retain 
the design order property. 




asymptotic 

range 

- 


° O o 

°o : 

o 

coarse 


1 . 6 1 „ , , 

10 ” 10 10 ” 10 ” 10 


1/N 


IV. 2-D Euler Equations 

A legacy CFD code used in this section. The code was a multi-block, 
finite volume code that solves the three-dimensional compressible Navier- 
Stokes equations. The inviscid terms are upwind-biased spatial differenced. 

Solutions are developed using the flux-difference splitting procedure of Roe 7 
with MUSCL 8 (Monotone Upstream-center Scheme for Conservative Laws) first order, second or higher 
order differencing. The solution is advanced in time with a three-factor approximately factored scheme. For 
the initial verification of the code the 2-D Euler equations were solved, reducing the governing equations to 
eqs. 33 through 37. 


Figure 8. 2-D Burger’s Equa- 

tion, Observed order vs. 1 /iV, 
double precision. 0> Code ; — 
0(Ai 2 ). 


dp dpu dpv 
dt + dx + dy ~ 

(33) 

lir + k“ 2 + p l + i ; [ p “”l = s - m 

(34) 

^ + ^M + §^ 2 +i<]=s,«. 

(35) 

d d d 

Ht ( PE ) + ~dx^ U ( PE +P ^ + Hy [ V ( PE + = S *nergy,MS 

(36) 

P = ( 7 - 0 (p e ^ \p (u 2 + w 2 )') 

(37) 


Forcing (or MS source) terms are generated from substituting the manufactured solutions in to the steady 
state forms of equations 33 through 36. 


Scant, 


nuity ,MS 


tyM’ 


d d 

Su.ms = - — [pUiUj + p5ij] , S energyMS = - — [uj ( pE + p)] (38) 


dx i 


dxj 


u 3 ^^3 

The specific form of the manufactured solutions for the Navier-Stokes equations using density, velocity, 
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and pressure are shown in Equations 39 through 42. 


P(x , y) = Po + Px sin + py cos (39) 

u(x,y) = u 0 + u x sin ^ au ^ 7rx ^ + Uy cos (40) 

v(x,y) = v 0 + v x cos ^ a ™ nx 'j + Vy s i n (41) 

p(x, y) = po + Px cos ^ a P^ nx ^j + Py s i n ( 4 2) 


^continuity , ms ; <S momerlt „ miMS , and 5 ener9!/iMS are calculated by expanding eqs. 38 through 38 with eqs. 39 through 
42 using the algebraic system Maple. c 


IV. A. Solution Process 


(a) |i?| vs. 
mesh. 


Solutions were initialized with the MS functions, 
though this is not a preferred mode of operation 
as discussed in Knupp & Salari. 1 The solution dis- 
cretization error is calculated each iteration for the 
primitive variables p, u, v and p after the solution 
vector is updated. 

The computational domain was two dimensional 
with an even spaced mesh in the vertical and hor- 
izontal directions. The domain was one unit in 
length in both dimensions with the reference lengths 
L x = L y = 1. Periodic boundary conditions were 
imposed on the four edges of the domain. Euler 
flow solutions were run with no flux limiting. 

The solution residual and discretization error 
with iteration number for the 10 x 10 mesh for 
a 0.1% magnitude manufactured solution, for ex- 
ample, Pol Px = Pol Py = 0.001 and all solutions 
were run with a CFL= 8. The results for the 
first-order discretization scheme are shown in fig- 
ures 9(a) and 9(b). For these parameters itera- 
tive convergence to machine zero was obtained, fig- 
ure 9(a). Similarly the discretization error remained 
unchanged for several thousand iterations, though 

L 2 „ converged in fewer iterations than L 2 p . The ap- ( c ) Lz ’P vs - h 

- ; o, p- ; — 



Iteration, 10 x 10 (b) L 2 
mesh ; 
P- 


Iteration, 10 x 10 
P 5 5 u 5 — 7 




; Code O, p 
O(Air). 


□ , (d) Order vs. h 
□ , u ; O, p. 


Code O, p 


Figure 9. First-order scheme verification, 0.1% magnitude MS, Euler flow, 
CFL=8, no limiting. 


propriate first order behavior was observed for these 
parameters as shown in figure 9(c) and figure 9(d). 

Manufactured solutions were repeated using a 
higher order option, k = 1/3, in the MUSCL 

scheme. Upon initial inspection of the results, figure 10(c), the code appeared to be verifying 3 rd order- 
accurate, though there was an unexpected shift in the L^norm values for h = 2.25 and h = 3.375 (and to 
a lesser degree h = 6.75). The u-component of velocity and pressure displayed similar trends as density 
comparing the □ and O with the O symbols in figure 10(c). 

The discretization errors with solution iteration of density, u- velocity and pressure for each mesh densities 
are plotted in figure 11. It appears that the deviation of the I^norm for the density and velocity is due to 
the non-convergence of the discretization error at the higher mesh densities. The discretization error of the 
pressure was the only apparently converged quantity when only the discretization error is inspected. 


c Maple is a symbol computation system. Further information can be found at http://www.maplesoft.com/ 
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(a) \R\ vs. Iteration, 10 x 10 (b) L 2 vs. Iteration, 10 x 10 

mesh mesh ; , p ; , u ; — • — , 

V 




(c) L/ 2 ,p vs. h ; Code O , p ; □ , (d) Order vs. h ; Code O, p ; 
u ; O, p ; , 0( Ax 3 ). □, u ; O, p. 


Figure 10. Higher-order MUSCL scheme ( k , = 1/3) verification, 0.1% 
magnitude MS, Euler flow, CFL=8, no limiting. 




(a) L 2 ,p vs. Iteration 


(b) h 2 ,u vs. Iteration 


(c) L/ 2 ,p vs. Iteration 


Figure 11. Convergence of discretization error with mesh density, 0.1% magnitude MS, Euler flow, CFL=8, 
k , = 1/3, no limiting. ; , 10 x 10 ; , 20 x 20 ; — 40 x 40 ; — 80 x 80 ; x, 120 x 120 

Only for the coarser mesh densities, ( ) 10 x 10 and ( ) 20 x 20, are iteratively converged in the 

solution residual plots, fig. 12. Potentially, after a sufficiently large number of iterations, the higher density 
meshes would be iteratively converged in the discretization error as well as the solution residual, but the 
wall-clock execution times would become unrealistic for routine testing. Table 9 lists the approximate run 
times for each of the mesh levels with MS. A per-iteration timing was roughly 20 milliseconds-per-point. 

There are several (at this time unresolved) possibilities that could prevent the method from verifying. 
Problems could exist with, but not limited to, the implementation of the manufactured solutions in the code, 
the construction of the higher-order fluxes, the construction of the implicit terms or boundary conditions. 


V. Summary 

Verification is an important step in the evaluation of a numerical method and the concept is fairly 
straightforward. The implementation and analysis though are often problematic and require understanding 
the method of verification as well as the method being verified. 

The method of manufactured solutions was applied to “textbook” and “real world” problems. Arithmetic 
precision was found to be a relevant factor for even a simple textbook problem. Additionally the design of 
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Table 9. Timing for Manufactured Solutions Using 2-D Euler Method 


Mesh 

Iter, time 
(1/sec) 

Iterations 

Total 

time 

10 x 10 

0.002 

5000 

10 sec. 

20 x 20 

0.008 

50000 

7 min. 

40 x 40 

0.032 

150000 

1.3 hr. 

80 x 80 

0.128 

150000 

5.3 hr. 

120 x 120 

0.288 

150000 

12 hr. 


the discretization template influenced the order of accuracy analysis as did the mesh form and form of the 
manufactured solutions, though the departure from design order occurred in specific detail not in general 
trends. 

With increasing size and complexity of the method to be verified, 
potentially the parameters that require analysis and/or inspection need 
to be increased as well. In the final example of the verification of the 
Euler equations in the Navier-Stokes solver, the non-convergence of the 
discretization error could have been over-looked if either only the behavior 
of the pressure was analyzed or the solutions stopped prematurely. Also, 
decisions must be made as to the time and resources to be devoted the 
verification of a code. The larger mesh densities required more resources 
and the run times quickly become unrealistic for efficient or routine re- 
verification processes. 

A. Three point discretization in transformed space 



An expression for the truncation error of a first derivative expression 
of / in the transformed coordinate £ was derived in Hoffman. 6 Using 
that methodology the truncation error of the second derivative of / in the 
transformed space is derived. Consider 


Figure 12. Variation of solution 
residual with mesh, MUSCL scheme 
k — 1/3, 0.1% magnitude MS, Euler 

flow, CFL=8, no limiting. ; , 10 

x 10 ; , 20 x 20 ; — • — •, 40 x 40 

80 x 80 ; x, 120 x 120 


/ = (j)(x) = </>{g(Q) = V’(C) (Al) 

where the physical, x, and transformed, £, coordinates are functions of 
the form, 


£ = h{x), x = g(0 

The second derivative f xx in the transformed coordinates is 

(A2) 

fxx = ('0(^))xrc 

(A3) 

= £x 

(A4) 

= 

(A5) 


with the metrics for this 1-D problem reducing to 


x i\i o ( Xi + 1 x i—l) ) L — 

n 2 n x j 




As defined previously Ax+ = Xi+\ — Xi and Ax- = Xi — Xi-i, so that 

Xi+i - Xi - 1 = ( Ax + - Ax _) 
resulting in the expression for the metric !; x - 

2 


L: — 


(Ax + — Ax-) 


(A6) 


(A7) 

(A8) 
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A discretized equation for f xx can be derived from either Eq. A4 or Eq. A5. Expansion of Eq. A4 about 
the half nodes is written as 


fxx |j — (OV’OfJ . 

= (^L+l/2^« 


i+1/2 ^ x li-l/2^ li-1/2; 

AC = A£ + = 6+1 - & = A£_ = & - 6-1 


(A9) 

(A10) 

(All) 


where the value of the metrics at the half nodes x i+1 /2 and x^_ i/ 2 are taken as the numerical average of the 
metrics at the neighboring nodes 


^ L— 1/2 — 2 (AL + ^L-0 

and the derivatives V’c| i + 1 / 2 are 


_ 1 

*+i/2 _ 2 


(&L+1 + ^ x li) 


V’e 


i+1/2 


lfi + 1 - ifi 


V'el 


- 1/2 


ifi - V’i-l 

A? 


(A12) 


(A13) 


Cell-centered finite volume computational methods typically define the metrics of the mesh at nodes between 
cell centers eliminating the need for Eq. A12. 

As an aside, if instead we start with equations A14, A15, and A16 in the transformed space and substitute 
in to Eq. A5, 


V»« = 

6e = 
Oj = 


V’i+i - 2 ^i + fa- 1 

AS 2 

if i+1/2 ~ if i— 1/2 

I i+1/2 ~~ li-1/2 




we can write 


fxx |j — (OjV’ff + OV^Odi 


t2 | ifi+l-Zlfi+lfi-l , ^ | if i+1/2 if i— 1/2 f | 

Sx|,; + Sx|,; 


A£ 2 


A£ 


Expanding Of and assuming 


if i+1/2 = ^ (lfi+ 1 + ifi) , if i— 1/2 = \{lfi + V’i-l) 


(A14) 

(A15) 

(A16) 

(A17) 

(A18) 

(A19) 


the final form of the second derivative in the transformed space is equivalent to the previous expression 
written for f xx though the truncation error differs between the two methods. 

The truncation error associated with Eq. A10 can be derived by expanding if in the transformed space 
and substituting in to Eq. A5 using the discrete forms in Eqs. A14 and A15. Performing the same Taylor 
series expansion, here on 6(0 1 


V’i+i = ifi ± 6f g Off ± gV’fff + ^6ffff 


AO 


± T ^tel ^ 5 + Ji " ; ± ° ( A2 °) 

Substituting Eq. A20 in to Eq. A18 averaging if using the Eq. A19 

fxx\ i = (06ff + 06f0f)* + + g06fff0f ] A f 2 + (gg g£*V’ffffff + Y^&’V’fffffOf] AO 


1 

5040 


(06fffffff0f)i AO + 0 ( AO) 


(A21) 
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So that the expression 


fxx ■ — 


Z \£ 2 L 


- 1pi)£ x | j+l/2 - {ipi ~ V’*-l)^x| i _ 1/2 


(A22) 


has the truncation error 
1 


fx 


ITE 


+ g ZxiPttt&ctJ -U' 2 + ) -4' 


+ 


1 


5040 


A£ 6 + O (A£ 8 ) (A23) 


The alternative derivation expands the Taylor series about the half nodes i ± 1/2, 

= * ±4 (f ) + >4 (f r ± (#)’■ + (# y ± ° < am > 

and substituting in to Eq. A18 (written only to 0 (A£ 4 ))- 

fxx Ij = (4V’« + £*44*0* + ) (A25) 

The difference between Eqs. A24 and A25 is the smaller numerical coefficient in the CxV’CttC.xg term which 
slightly shifts the level of the truncation error but does not alter the order of the term. 

The discrete equation for f xx in the transformed space is formally second order in Z\£ for either method 
of forming if) x . 


B. 1-D Mesh Parameters 
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Table 10. Geometric Series Mesh Parameters 


Grid number 

N 

O.OOOIAejv 

r 

Grid number 

N 

O.OOOIAejv 

r 

1 

5 

1.667 

xlO -5 

8.8139 

16 

1851 

5.400 

xl0~ 8 

1.0063 

2 

7 

1.250 

xl0~ 5 

4.8543 

17 

2776 

3.601 

xlO -8 

1.0042 

3 

10 

9.091 

xlO -6 

3.0691 

18 

4164 

2.401 

xlO -8 

1.0028 

4 

15 

6.250 

xl0~ 6 

2.1311 

19 

6246 

1.601 

xlO -8 

1.0019 

5 

22 

4.348 

xlO" 6 

1.6823 

20 

9369 

1.067 

xlO" 8 

1.0012 

6 

33 

2.941 

xlO" 6 

1.4175 

21 

14053 

7.115 

xlO" 9 

1.0008 

7 

49 

2.000 

xlO" 6 

1.2661 

22 

21079 

4.743 

xlO -9 

1.0006 

8 

73 

1.351 

xlO" 6 

1.1722 

23 

31618 

3.163 

xlO" 9 

1.0004 

9 

109 

9.091 

xlO" 7 

1.1125 

24 

47427 

2.109 

xlO" 9 

1.0002 

10 

163 

6.098 

xlO" 7 

1.0740 

25 

71140 

1.406 

xlO" 9 

1.0002 

11 

244 

4.082 

xlO" 7 

1.0489 

26 

106710 

9.371 

xlO" 10 

1.0001 

12 

366 

2.725 

xlO" 7 

1.0324 

27 

160065 

6.247 

xlO" 10 

1.0001 

13 

549 

1.818 

xlO" 7 

1.0215 

28 

240097 

4.165 

xlO" 10 

<1.0001 

14 

823 

1.214 

xlO" 7 

1.0143 

29 

360145 

2.777 

xlO -10 

<1.0001 

15 

1234 

8.097 

xlO" 8 

1.0095 

30 

540217 

1.851 

xlO -10 

<1.0001 


Table 11. Elliptic DE Mesh Parameters 


Grid number 

N 

Axn 

Ax o 

Grid number 

N 

Axn 

Ax o 

1 

5 

1.6667E-01 

1.7510E-05 

16 

1851 

5.3996E-04 

3.6474E-07 

2 

7 

1.2500E-01 

4.9805E-05 

17 

2776 

3.6010E-04 

2.4303E-07 

3 

10 

9.0909E-02 

5.5979E-05 

18 

4164 

2.4010E-04 

1.6195E-07 

4 

15 

6.2500E-02 

4.5053E-05 

19 

6246 

1.6008E-04 

1.0793E-07 

5 

22 

4.3478E-02 

3.2401E-05 

20 

9369 

1.0672E-04 

7.1938E-08 

6 

33 

2.9412E-02 

2.1768E-05 

21 

14053 

7.1154E-05 

4.7953E-08 

7 

49 

2.0000E-02 

1.4530E-05 

22 

21079 

4.7438E-05 

3.1967E-08 

8 

73 

1.3514E-02 

9.6318E-06 

23 

31618 

3.1627E-05 

2.1310E-08 

9 

109 

9.0909E-03 

6.3768E-06 

24 

47427 

2.1085E-05 

1.4206E-08 

10 

163 

6.0976E-03 

4.2254E-06 

25 

71140 

1.4057E-05 

9.4705E-09 

11 

244 

4.0816E-03 

2.8036E-06 

26 

106710 

9.3711E-06 

6.3136E-09 

12 

366 

2.7248E-03 

1.8600E-06 

27 

160065 

6.2474E-06 

4.2090E-09 

13 

549 

1.8182E-03 

1.2358E-06 

28 

240097 

4.1650E-06 

2.8060E-09 

14 

823 

1.2136E-03 

8.2249E-07 

29 

360145 

2.7767E-06 

1.8706E-09 

15 

1234 

8.0972E-04 

5.4769E-07 

30 

540217 

1.8511E-06 

1.2471E-09 
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